// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ***************************************************************************
//                            avtOpenPMDFileFormat.C
// ***************************************************************************
#ifdef _WIN32
#include <io.h>
#else
#include <unistd.h>
#endif

#include <avtOpenPMDFileFormat.h>

#include <vector>
#include <stdlib.h>
#include <stdio.h>
#include <stdio.h>
#include <iostream>
#include <string>

#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPoints.h>
#include <vtkVertex.h>

#include <avtDatabaseMetaData.h>

#include <Expression.h>
#include <DebugStream.h>

#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>
#include <InvalidFilesException.h>
#include <InvalidTimeStepException.h>
#include <InvalidDBTypeException.h>

using     std::string;


// ****************************************************************************
//  Method: avtOpenPMDFileFormat constructor
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************

avtOpenPMDFileFormat::avtOpenPMDFileFormat(const char *filename)
    : avtMTMDFileFormat(filename)
{
    debug5 << "avtOpenPMDFileFormat::avtOpenPMDFileFormat "
           << filename
           << endl;

    // Boolean to check that the initialization has been done
    this->initialized = false;
    // By default parallel is off
    this->parallel = false;
}


// ****************************************************************************
//  Method: ~avtOpenPMDFileFormat Destructor
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************
avtOpenPMDFileFormat::~avtOpenPMDFileFormat()
{
    debug5 << "avtOpenPMDFileFormat::~avtOpenPMDFileFormat" << endl;
}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::Initialize
//
//  Initialization: the openPMD file is read and analyze
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************
void
avtOpenPMDFileFormat::Initialize()
{
    debug5 << "avtOpenPMDFileFormat::Initialize" << endl;
    if(!this->initialized)
    {
        H5Eset_auto(H5E_DEFAULT, NULL, NULL);

        // Open the OpenPMD file
        openPMDFile.OpenFile(filename);

        // Get the root attributes
        openPMDFile.ScanFileAttributes();

        // Scan the iterations
        openPMDFile.ScanIterations();

        // For each iteration, scan the fields
        openPMDFile.ScanFields();

        // For each iteration, scan the particle groups
        openPMDFile.ScanParticles();

        // Number of MPI tasks
        this->numTasks = PAR_Size();
        // If more than one task, we use a parallel treatment
        if (this->numTasks > 1)
        {
            this->parallel = true;
            debug5 << " Parallel with " << this->numTasks << " tasks"<< endl;
        }

        // Initialization has been done
        this->initialized = true;
    }
}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::GetNTimesteps
//
//  Purpose:
//      Tells the rest of the code how many timesteps there are in this file.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************

int
avtOpenPMDFileFormat::GetNTimesteps(void)
{
    debug5 << "avtOpenPMDFileFormat::GetNTimesteps" << endl;

    // If not initialized, initialize...
    if(!this->initialized)
        this->Initialize();

    // Get the number of time steps read in the OpenPMD file
    int NTimesteps = this->openPMDFile.GetNumberIterations();

    return NTimesteps;
    //return 1;
}

// ****************************************************************************
//  Method: avtOpenPMDFileFormat::GetCycles
//
//  Purpose:
//      Get all cycles in the openPMD file.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************
void
avtOpenPMDFileFormat::GetCycles(std::vector<int> &cycles)
{

    int cycle;

    // Shortcut pointer to the vector iterations
    vector<PMDIteration> * iterations = &(openPMDFile.iterations);

    // Loop over the iterations
    for (std::vector<PMDIteration>::iterator it = iterations->begin() ;
         it != iterations->end(); ++it)
    {
        cycle = atoi(it->name.c_str());

        // Store the cycles in the vector.
        cycles.push_back(cycle);
    }

}

// ****************************************************************************
//  Method: avtOpenPMDFileFormat::GetTime
//
//  Purpose:
//      Get all times in the openPMD file.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************
void
avtOpenPMDFileFormat::GetTimes(std::vector<double> &times)
{
    debug5 << "avtOpenPMDFileFormat::GetTimes" << endl;
    // Shortcut pointer to the vector iterations
    vector<PMDIteration> * iterations = &(openPMDFile.iterations);

    // Loop over the iterations
    for (std::vector<PMDIteration>::iterator it = iterations->begin() ;
         it != iterations->end(); ++it)
    {
        // Store the times in the vector.
        times.push_back(it->time*it->timeUnitSI);
    }
}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
// ****************************************************************************

void
avtOpenPMDFileFormat::FreeUpResources(void)
{
    debug5 << "avtOpenPMDFileFormat::FreeUpResources" << endl;
}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
//  Modifications:
//  - Nov 25 2016 - Mathieu Lobet - Particle can now be read by block
//                                  depending of the number of tasks.
//
// ****************************************************************************

void
avtOpenPMDFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md,
                                               int timeState)
{
    debug5 << "avtOpenPMDFileFormat::PopulateDatabaseMetaData(timeState="
         << timeState << ")"
         << endl;

    int i;
    string buffer;
    string bufferMeshName;

    // ________________________________________________________
    // FIELDS

    // Creation of the meshes and the scalars from the vector of field objects

    // Shortcut pointer to the fields
    vector<PMDField> * fields = &(openPMDFile.iterations[timeState].fields);

    // Iteration over the fields
    for (std::vector<PMDField>::iterator field = fields->begin() ;
         field != fields->end(); ++field)
    {

        // We first create a new mesh for the current field
        // Each field has its own mesh
        avtMeshMetaData *mmd = new avtMeshMetaData;
        // Create the mesh name
        buffer = "Fields/";
        buffer += field->name;
        buffer += "_mesh";
        mmd->name = buffer.c_str();
        // Spatial dimension
        mmd->spatialDimension = field->ndims;
        // Topological Dimension
        mmd->topologicalDimension = field->ndims;
        // Mesh type
        // Cartesian
        if (field->geometry == "cartesian")
        {
            mmd->meshType = AVT_RECTILINEAR_MESH;
        }
        // Curvilinear, thetaMode
        else if (field->geometry == "thetaMode")
        {
            mmd->meshType = AVT_CURVILINEAR_MESH;
        }

        // Labels

        if (field->dataOrder == "C")
        {
            mmd->xLabel = field->axisLabels[2];
            mmd->yLabel = field->axisLabels[1];
            mmd->zLabel = field->axisLabels[0];
        }
        else if (field->dataOrder == "Fortran")
        {
            mmd->xLabel = field->axisLabels[0];
            mmd->yLabel = field->axisLabels[1];
            mmd->zLabel = field->axisLabels[2];
        }

        // Number of blocks
        mmd->numBlocks = this->numTasks;
        // Add mesh
        md->Add(mmd);

        // Then we create a scalar for the current field
        avtScalarMetaData *smd = new avtScalarMetaData;
        // Create the variable name
        buffer = "Fields/";
        buffer += field->name;
        smd->name = buffer.c_str();
        // And select the right mesh
        buffer = "Fields/";
        buffer += field->name;
        buffer += "_mesh";
        smd->meshName = buffer.c_str();
        // Node or cell centered
        smd->centering = AVT_NODECENT;
        // Units
        smd->hasUnits = true;
        smd->units = field->unitsLabel;
        // Add the scalars
        md->Add(smd);

    }

    // ____________________________________________________________
    // FIELD GROUP

    // Shortcut pointer to the field groups
    vector<fieldGroupStruct> * fieldGroups
                             = &(openPMDFile.iterations[timeState].fieldGroups);

    // Iteration over the field group for additional components
    // computed with multiple datasets
    for (std::vector<fieldGroupStruct>::iterator
         fieldGroup = fieldGroups->begin();
         fieldGroup != fieldGroups->end();
         ++fieldGroup)
    {

      // Additional field components for the theta mode
      if (fieldGroup->geometry == "thetaMode")
      {
          // check if components r and theta are available
          // If True, we compute x and y from r and theta.
          if ((fieldGroup->thetaComponents[0]>-1) &&
              (fieldGroup->thetaComponents[1]>-1))
          {

            // We create a scalar for the x components
            // that will be computed from r and theta
            avtScalarMetaData *smdEx = new avtScalarMetaData;
            // Create the variable name
            buffer = "Fields/";
            buffer += fieldGroup->name;
            buffer += "/x";
            smdEx->name = buffer.c_str();
            // And select the right mesh
            // We use the same mesh as for the r component
            buffer = "Fields/";
            buffer += fieldGroup->name;
            buffer += "/r_mesh";
            smdEx->meshName = buffer.c_str();
            // Node or cell centered
            smdEx->centering = AVT_NODECENT;
            // Units
            smdEx->hasUnits = true;
            // We use the unit label of the r component
            smdEx->units = openPMDFile.iterations[timeState].
                         fields[fieldGroup->thetaComponents[0]].unitsLabel;
            // Add the scalars
            md->Add(smdEx);

            // We create a scalar for the x components
            // that will be computed from r and theta
            avtScalarMetaData *smdEy = new avtScalarMetaData;
            // Create the variable name
            buffer = "Fields/";
            buffer += fieldGroup->name;
            buffer += "/y";
            smdEy->name = buffer.c_str();
            // And select the right mesh
            // We use the same mesh as for the r component
            buffer = "Fields/";
            buffer += fieldGroup->name;
            buffer += "/r_mesh";
            smdEy->meshName = buffer.c_str();
            // Node or cell centered
            smdEy->centering = AVT_NODECENT;
            // Units
            smdEy->hasUnits = true;
            // We use the unit label of the r component
            smdEy->units = openPMDFile.iterations[timeState].
                         fields[fieldGroup->thetaComponents[0]].unitsLabel;
            // Add the scalars
            md->Add(smdEy);

          }
      }
    }


    // ________________________________________________________
    // PARTICLES

    // Shortcut pointer to the particles
    vector<PMDParticle> * particles =
                            &(openPMDFile.iterations[timeState].particles);

    // Iteration over the particles
    for (std::vector<PMDParticle>::iterator particle = particles->begin() ;
         particle != particles->end(); ++particle)
    {

        // Add a point mesh to the metadata for the particles
        avtMeshMetaData *mmd = new avtMeshMetaData;
        // Mesh name
        bufferMeshName = "Particles/";
        bufferMeshName += particle->name;
        bufferMeshName += "/position_mesh";
        mmd->name = bufferMeshName.c_str();
        // Spatial and topological dimension
        mmd->spatialDimension = 3;
        mmd->topologicalDimension = 0;
        // Type of mesh, here point mesh
        mmd->meshType = AVT_POINT_MESH;
        // The number of blocks is equal to the number of parallel tasks
        mmd->numBlocks = this->numTasks;
        md->Add(mmd);

        // We create a scalar metadata for each scalar in scalarDataSets
        for (i=0;i<particle->GetNumScalarDatasets();i++)
        {

            // Then we create a scalar for the current field
            avtScalarMetaData *smd = new avtScalarMetaData;
            // Create the scalar name
            buffer = "Particles/";
            buffer += particle->scalarDataSets[i].name;
            smd->name = buffer.c_str();

            // And select the right mesh
            smd->meshName = bufferMeshName.c_str();
            // Node or cell centered
            smd->centering = AVT_NODECENT;

            // Units
            smd->hasUnits = true;
            smd->units = particle->scalarDataSets[i].unitLabel;
            // Add the scalars
            md->Add(smd);

        }

        // We create a vector metadata for each scalar in scalarDataSets
        for (i=0;i<particle->GetNumVectorDatasets();i++)
        {
            avtVectorMetaData *vmd = new avtVectorMetaData;
            // Create the scalar name
            buffer = "Particles/";
            buffer += particle->vectorDataSets[i].name;
            vmd->name = buffer.c_str();
            // And select the right mesh
            vmd->meshName = bufferMeshName.c_str();
            // Node or cell centered
            vmd->centering = AVT_ZONECENT;
            // Units
            vmd->hasUnits = true;
            vmd->units = particle->vectorDataSets[i].unitLabel;
            // Vector dimension
            vmd->varDim = 3;
            md->Add(vmd);
        }

        if (particle->momentumAvailable)
        {
            // Particle Expressions
            Expression *e1 = new Expression;
            // Create the name
            buffer = "Particles/";
            buffer += particle->name;
            buffer += "/gamma";
            e1->SetName(buffer);

            // Creation of the expression definition
            char definition[256];

            // If we have 3 momentum components
            if (particle->numDimsMomenta==3)
            {
                // Create the definition
                char px[128], py[128], pz[128];

                // Expression
                sprintf(px,"Particles/%s/momentum/x",particle->name);
                sprintf(py,"Particles/%s/momentum/y",particle->name);
                sprintf(pz,"Particles/%s/momentum/z",particle->name);

                sprintf(definition,
                        "sqrt(1.0+(<%s>^2+<%s>^2+<%s>^2)/(299792458*%e)^2)",
                        px,py,pz,particle->mass);
            }
            // If we have only two momentum components, we suppose
            // that the dimension is 2
            else if (particle->numDimsMomenta==2)
            {
                // Create the definition
                char px[128], pz[128];

                // Expression: this is only valid if z is the main direction
                sprintf(px,"Particles/%s/momentum/x",particle->name);
                sprintf(pz,"Particles/%s/momentum/z",particle->name);

                sprintf(definition,
                        "sqrt(1.0+(<%s>^2+<%s>^2)/(299792458*%e)^2)",
                        px,pz,particle->mass);
            }
            // Invalid momentum dimension
            else
            {
                debug5 << " The gamma expression can not be created for "
                        "such a momentum dimension: "
                       << particle->numDimsMomenta << endl;
            }

            e1->SetDefinition(definition);
            // Final type
            e1->SetType(Expression::ScalarMeshVar);
            // Units
            //e1->hasUnits = true;
            //e1->units = "";
            // Expression is visible to users
            e1->SetHidden(false);
            // Add expression to the list of metadata
            md->AddExpression(e1);
        }
    }




    //
    // CODE TO ADD A MESH
    //
    // string meshname = ...
    //
    // AVT_RECTILINEAR_MESH, AVT_CURVILINEAR_MESH, AVT_UNSTRUCTURED_MESH,
    // AVT_POINT_MESH, AVT_SURFACE_MESH, AVT_UNKNOWN_MESH
    // avtMeshType mt = AVT_RECTILINEAR_MESH;
    //
    // int nblocks = YOU_MUST_DECIDE;
    // int block_origin = 0;
    // int spatial_dimension = 2;
    // int topological_dimension = 2;
    // double *extents = NULL;
    //
    // Here's the call that tells the meta-data object that we have a mesh:
    //
    // AddMeshToMetaData(md, meshname, mt, extents, nblocks, block_origin,
    //                   spatial_dimension, topological_dimension);
    //

    //
    // CODE TO ADD A SCALAR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddScalarVarToMetaData(md, varname, mesh_for_this_var, cent);
    //

    //
    // CODE TO ADD A VECTOR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    // int vector_dim = 2;
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddVectorVarToMetaData(md, varname, mesh_for_this_var, cent,vector_dim);
    //

    //
    // CODE TO ADD A TENSOR VARIABLE
    //
    // string mesh_for_this_var = meshname; // ??? -- could be multiple meshes
    // string varname = ...
    // int tensor_dim = 9;
    //
    // AVT_NODECENT, AVT_ZONECENT, AVT_UNKNOWN_CENT
    // avtCentering cent = AVT_NODECENT;
    //
    //
    // Here's the call that tells the meta-data object that we have a var:
    //
    // AddTensorVarToMetaData(md, varname, mesh_for_this_var, cent,tensor_dim);
    //

    //
    // CODE TO ADD A MATERIAL
    //
    // string mesh_for_mat = meshname; // ??? -- could be multiple meshes
    // string matname = ...
    // int nmats = ...;
    // vector<string> mnames;
    // for (int i = 0 ; i < nmats ; i++)
    // {
    //     char str[32];
    //     sprintf(str, "mat%d", i);
    //     -- or --
    //     strcpy(str, "Aluminum");
    //     mnames.push_back(str);
    // }
    //
    // Here's the call that tells the meta-data object that we have a mat:
    //
    // AddMaterialToMetaData(md, matname, mesh_for_mat, nmats, mnames);
    //
    //
    // Here's the way to add expressions:
    //Expression momentum_expr;
    //momentum_expr.SetName("momentum");
    //momentum_expr.SetDefinition("{u, v}");
    //momentum_expr.SetType(Expression::VectorMeshVar);
    //md->AddExpression(&momentum_expr);
    //Expression KineticEnergy_expr;
    //KineticEnergy_expr.SetName("KineticEnergy");
    //KineticEnergy_expr.SetDefinition("0.5*(momentum*momentum)/(rho*rho)");
    //KineticEnergy_expr.SetType(Expression::ScalarMeshVar);
    //md->AddExpression(&KineticEnergy_expr);
    //
}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      timestate   The index of the timestate.  If GetNTimesteps returned
//                  'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
//  Modifications:
//
//  - Nov 14 2016 - Mathieu Lobet - Add parallel treatment of the fields
//
//  - Nov 25 2016 - Mathieu Lobet - Add parallel treatment of the particles
//
//  - Feb  8 2016 - Mathieu Lobet - Add double precision for particles
//                                  and fields
//
// ****************************************************************************

vtkDataSet *
avtOpenPMDFileFormat::GetMesh(int timestate, int domain, const char *meshname)
{

    debug5 << "avtOpenPMDFileFormat::GetMesh(timestate="
         << timestate
         << ",domain=" << domain
         << ")" << endl;

    int             ndims;
    int             dims[3];
    int             i,j,k,l;
    int             i0,i2;
    int             err;
    int             id;
    int             nnodes;
    int             nmodes;
    int             dataSize;
    string          bufferMeshName;
    char            bufferDataSetName[128];
    float           factor;
    float           r,theta,dtheta;
    vtkFloatArray*  coords[3];

    // ________________________________________________________
    // FIELDS

    // We look for the corresponding mesh of name meshname
    // in the vector of fields

    // Shortcut pointer to the fields
    vector<PMDField> * fields = &(openPMDFile.iterations[timestate].fields);

    // Iteration over the fields
    for (std::vector<PMDField>::iterator field = fields->begin() ;
         field != fields->end(); ++field)
    {

        // Current field meshname
        bufferMeshName = "Fields/";
        bufferMeshName += field->name;
        bufferMeshName += "_mesh";

        // We have found the corresponding mesh in the vector of fields
        if (strcmp(meshname, bufferMeshName.c_str()) == 0)
        {

            if (field->geometry == "cartesian")
            {
                // Read the ndims and number of X,Y,Z nodes from file.
                ndims = field->ndims;

                // 3D
                if (ndims==3)
                {

                    // Data order
                    if (field->dataOrder=="C")
                    {
                        i0 = 2;
                        i2 = 0;
                    }
                    else if (field->dataOrder == "Fortran")
                    {
                        i0 = 0;
                        i2 = 2;
                    }

                    // Treatment of the file in parallel
                    if (this->parallel)
                    {

                        // Structure to store the block properties
                        fieldBlockStruct fieldBlock;

                        // We get the block properties
                        err = field->GetBlockProperties(this->numTasks,
                                                        domain,
                                                        &fieldBlock);

                        debug5    << "fieldBlock.minNode[0]: "
                                << fieldBlock.minNode[0]
                                << " " << fieldBlock.nbNodes[0]
                                << " " << fieldBlock.maxNode[0] << endl;

                        debug5    << "fieldBlock.minNode[1]: "
                                << fieldBlock.minNode[1]
                                << " " << fieldBlock.nbNodes[1]
                                << " " << fieldBlock.maxNode[1] << endl;

                        debug5    << "fieldBlock.minNode[2]: "
                                << fieldBlock.minNode[2]
                                << " " << fieldBlock.nbNodes[2]
                                << " " << fieldBlock.maxNode[2] << endl;

                        // Dimensions
                        dims[0] = fieldBlock.nbNodes[i0];
                        dims[1] = fieldBlock.nbNodes[1];
                        dims[2] = fieldBlock.nbNodes[i2];

                        // Read the X coordinates from the file.
                        coords[0] = vtkFloatArray::New();
                        coords[0]->SetNumberOfTuples(dims[0]);
                        float *xarray = (float *)coords[0]->GetVoidPointer(0);
                        for(i=fieldBlock.minNode[i0];
                            i<=fieldBlock.maxNode[i0];i++)
                        {
                            xarray[i - fieldBlock.minNode[i0]] =
                            ((i+field->gridPosition[i0])*field->gridSpacing[i0]
                            + field->gridGlobalOffset[i0])*field->gridUnitSI;
                        }

                        // Read the Y coordinates from the file.
                        coords[1] = vtkFloatArray::New();
                        coords[1]->SetNumberOfTuples(dims[1]);
                        float *yarray = (float *)coords[1]->GetVoidPointer(0);
                        for(i=fieldBlock.minNode[1];
                            i<=fieldBlock.maxNode[1];i++)
                        {
                            yarray[i - fieldBlock.minNode[1]] =
                            ((i+field->gridPosition[1])*field->gridSpacing[1]
                            + field->gridGlobalOffset[1])*field->gridUnitSI;
                        }

                        // Read the Z coordinates from the file.
                        coords[2] = vtkFloatArray::New();
                        coords[2]->SetNumberOfTuples(dims[2]);
                        float *zarray = (float *)coords[2]->GetVoidPointer(0);
                        for(i=fieldBlock.minNode[i2];
                            i<=fieldBlock.maxNode[i2];i++)
                        {
                            zarray[i - fieldBlock.minNode[i2]] =
                            ((i+field->gridPosition[i2])*field->gridSpacing[i2]
                            + field->gridGlobalOffset[i2])*field->gridUnitSI;
                        }

                    }
                    // Only one processor
                    else
                    {

                        dims[0] = field->nbNodes[i0];
                        dims[1] = field->nbNodes[1];
                        dims[2] = field->nbNodes[i2];

                        // Read the X coordinates from the file.
                        coords[0] = vtkFloatArray::New();
                        coords[0]->SetNumberOfTuples(dims[0]);
                        float *xarray = (float *)coords[0]->GetVoidPointer(0);
                        for(i=0;i<dims[0];i++)
                        {
                            xarray[i] =
                            ((i+field->gridPosition[i0])*field->gridSpacing[i0]
                            + field->gridGlobalOffset[i0])*field->gridUnitSI;
                        }

                        // Read the Y coordinates from the file.
                        coords[1] = vtkFloatArray::New();
                        coords[1]->SetNumberOfTuples(dims[1]);
                        float *yarray = (float *)coords[1]->GetVoidPointer(0);
                        for(i=0;i<dims[1];i++)
                        {
                            yarray[i] =
                            ((i+field->gridPosition[1])*field->gridSpacing[1]
                            + field->gridGlobalOffset[1])*field->gridUnitSI;
                        }

                        // Read the Z coordinates from the file.
                        coords[2] = vtkFloatArray::New();
                        coords[2]->SetNumberOfTuples(dims[2]);
                        float *zarray = (float *)coords[2]->GetVoidPointer(0);
                        for(i=0;i<dims[2];i++)
                        {
                            zarray[i] =
                            ((i+field->gridPosition[i2])*field->gridSpacing[i2]
                            + field->gridGlobalOffset[i2])*field->gridUnitSI;
                        }
                    }

                }
                // 2D
                else if (ndims==2)
                {

                    // Treatment of the file in parallel
                    if (this->parallel)
                    {

                        // Structure to store the block properties
                        fieldBlockStruct fieldBlock;

                        // We get the block properties
                        err = field->GetBlockProperties(this->numTasks,
                                                        domain,
                                                        &fieldBlock);

                        debug5    << "fieldBlock.minNode[0]: "
                                  << fieldBlock.minNode[0]
                                  << " " << fieldBlock.nbNodes[0]
                                  << " " << fieldBlock.maxNode[0] << endl;

                        debug5    << "fieldBlock.minNode[1]: "
                                  << fieldBlock.minNode[1]
                                  << " " << fieldBlock.nbNodes[1]
                                  << " " << fieldBlock.maxNode[1] << endl;

                        dims[0] = fieldBlock.nbNodes[0];
                        dims[1] = fieldBlock.nbNodes[1];
                        dims[2] = 1;

                        // Read the X coordinates from the file.
                        coords[0] = vtkFloatArray::New();
                        coords[0]->SetNumberOfTuples(dims[0]);
                        float *xarray = (float *)coords[0]->GetVoidPointer(0);
                        for(i=fieldBlock.minNode[0];
                            i<=fieldBlock.maxNode[0];i++)
                        {
                            xarray[i - fieldBlock.minNode[0]] =
                            ((i+field->gridPosition[0])*field->gridSpacing[0]
                            + field->gridGlobalOffset[0])*field->gridUnitSI;
                        }

                        // Read the Y coordinates from the file.
                        coords[1] = vtkFloatArray::New();
                        coords[1]->SetNumberOfTuples(dims[1]);
                        float *yarray = (float *)coords[1]->GetVoidPointer(0);
                        for(i=fieldBlock.minNode[1];
                            i<=fieldBlock.maxNode[1];i++)
                        {
                            yarray[i - fieldBlock.minNode[1]] =
                            ((i+field->gridPosition[1])*field->gridSpacing[1]
                            + field->gridGlobalOffset[1])*field->gridUnitSI;
                        }

                        // No Z coordinates
                        coords[2] = vtkFloatArray::New();
                        coords[2]->SetNumberOfTuples(1);
                        coords[2]->SetComponent(0, 0, 0.);

                    }
                    // Only one processor
                    else
                    {

                        dims[0] = field->nbNodes[0];
                        dims[1] = field->nbNodes[1];
                        dims[2] = 1;

                        // Read the X coordinates from the file.
                        coords[0] = vtkFloatArray::New();
                        coords[0]->SetNumberOfTuples(field->nbNodes[0]);
                        float *xarray = (float *)coords[0]->GetVoidPointer(0);
                        for(i=0;i<field->nbNodes[0];i++)
                        {
                            xarray[i] =
                            ((i+field->gridPosition[0])*field->gridSpacing[0]
                            + field->gridGlobalOffset[0])*field->gridUnitSI;
                        }

                        // Read the Y coordinates from the file.
                        coords[1] = vtkFloatArray::New();
                        coords[1]->SetNumberOfTuples(field->nbNodes[1]);
                        float *yarray = (float *)coords[1]->GetVoidPointer(0);
                        for(i=0;i<field->nbNodes[1];i++)
                        {
                            yarray[i] =
                            ((i+field->gridPosition[1])*field->gridSpacing[1]
                            + field->gridGlobalOffset[1])*field->gridUnitSI;
                        }

                        // No Z coordinates
                        coords[2] = vtkFloatArray::New();
                        coords[2]->SetNumberOfTuples(1);
                        coords[2]->SetComponent(0, 0, 0.);
                    }
                }
                // No recognized dimension
                else
                {
                    debug5 << " Dimension of this field "
                    "dataset is not recognized: " << ndims << endl;
                }

                // Create the vtkRectilinearGrid object and set its dimensions
                // and coordinates.
                vtkRectilinearGrid *grid = vtkRectilinearGrid::New();
                grid->SetDimensions(dims);
                grid->SetXCoordinates(coords[0]);
                coords[0]->Delete();
                grid->SetYCoordinates(coords[1]);
                coords[1]->Delete();
                grid->SetZCoordinates(coords[2]);
                coords[2]->Delete();

                // We return the generated grid
                return grid;
            }
            // Axicylindrical case
            else if (field->geometry == "thetaMode")
            {

                if (field->dataOrder == "C")
                {
                    // Read the ndims and number of x,r,theta nodes from file.
                    ndims = field->ndims;
                    dims[0] = field->nbNodes[2]; // x direction
                    dims[1] = field->nbNodes[1]; // r direction
                    dims[2] = field->thetaNbNodes; // Theta direction

                }
                else if (field->dataOrder == "Fortran")
                {
                    // Read the ndims and number of z,r,theta nodes from file.
                    ndims = field->ndims;
                    dims[0] = field->nbNodes[1]; // z direction
                    dims[1] = field->nbNodes[2]; // r direction
                    dims[2] = field->thetaNbNodes; // Theta direction
                }

                // Total number of nodes
                nnodes = dims[0] * dims[1] * dims[2];
                // Number of modes
                nmodes = field->nbModes;

                // Allocate coordinate arrays
                float *xarray = new float[nnodes];
                float *yarray = new float[nnodes];
                float *zarray = new float[nnodes];

                // Step in for the theta discretization
                dtheta = 2.*3.14159265359/(dims[2]-1);

                // Build xarray,yarray,zarray
                for(k = 0; k < dims[2]; ++k)
                for(j = 0; j < dims[1]; ++j)
                for(i = 0; i < dims[0]; ++i)
                {
                    // Absolute index
                    l = i + j*dims[0] + k*dims[0]*dims[1];
                    // Theta angle
                    theta = k*dtheta;
                    // Radius
                    r = (j*field->gridSpacing[0]
                      + field->gridGlobalOffset[0])*field->unitSI;
                    // Position of the point
                    xarray[l] = (i*field->gridSpacing[1]
                              + field->gridGlobalOffset[1])*field->gridUnitSI;
                    yarray[l] = (r*cos(theta));
                    zarray[l] = (r*sin(theta));
                }

                // Create the vtkStructuredGrid and vtkPoints objects.
                vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
                vtkPoints *points = vtkPoints::New();
                sgrid->SetPoints(points);
                sgrid->SetDimensions(dims);
                points->Delete();
                points->SetNumberOfPoints(nnodes);

                // Copy the coordinate values into the vtkPoints object.
                float *pts = (float *) points->GetVoidPointer(0);
                float *xc = xarray;
                float *yc = yarray;
                float *zc = zarray;
                if(ndims == 3)
                {
                    for(k = 0; k < dims[2]; ++k)
                    for(j = 0; j < dims[1]; ++j)
                    for(i = 0; i < dims[0]; ++i)
                    {
                        *pts++ = *xc++;
                        *pts++ = *yc++;
                        *pts++ = *zc++;
                    }
                }
                else if(ndims == 2)
                {
                    for(j = 0; j < dims[1]; ++j)
                    for(i = 0; i < dims[0]; ++i)
                    {
                        *pts++ = *xc++;
                        *pts++ = *yc++;
                        *pts++ = 0.;
                    }
                }
                // Delete temporary arrays.
                delete [] xarray;
                delete [] yarray;
                delete [] zarray;
                return sgrid;
            }
        }
    }

    // ________________________________________________________
    // PARTICLES

    // Shortcut pointer to the particles
    vector<PMDParticle> * particles =
    &(openPMDFile.iterations[timestate].particles);

    // Iteration over the particles
    for (std::vector<PMDParticle>::iterator particle = particles->begin() ;
         particle != particles->end(); ++particle)
    {
        // Mesh name
        bufferMeshName = "Particles/";
        bufferMeshName += particle->name;
        bufferMeshName += "/position_mesh";

        if (strcmp(meshname, bufferMeshName.c_str()) == 0)
        {

            // Read the dimension ndims
            ndims = particle->numDimsPositions;

            // Sanity check
            if ((ndims <= 1)||(ndims > 3))
            {
                debug5 << "The number of position datasets is invalid: "
                       << ndims << endl;
            }

            /* ____ Treatment of the file in parallel _____________________ */
            if (this->parallel)
            {

                // Creation of the block property object
                particleBlockStruct particleBlock;

                /* ____ Read the X coordinates from the file ______________ */

                // Dataset Id
                id = particle->positionsId[2];

                // Get the block properties for this dataSet
                particle->GetBlockProperties(id,
                                             this->numTasks,
                                             domain,
                                             &particleBlock);

                // Dataset path
                strcpy(bufferDataSetName,particleBlock.dataSetPath);

                // Multiplication factor
                factor = particle->scalarDataSets[id].unitSI;

                // DataSize
                dataSize = particle->scalarDataSets[id].dataSize;

                // Depending on the data size of this first dataSet,
                // we suppose that all data have the same type.

                // Simple precision float
                if (dataSize == 4)
                {
                    // Create the X coordinates array
                    float *xarray = new float[particleBlock.numParticles];

                    // Read the dataset
                    openPMDFile.ReadParticleScalarBlock(xarray,
                                                        &factor,
                                                        H5T_FLOAT,
                                                         &particleBlock);

                    /* ____ Read the Y coordinates from the file _________ */

                    // Dataset Id
                    id = particle->positionsId[1];

                    // Get the block properties for this dataSet
                    particle->GetBlockProperties(id, this->numTasks,
                                                 domain,
                                                 &particleBlock);

                    // Dataset path
                    strcpy(bufferDataSetName,particleBlock.dataSetPath);

                    // Multiplication factor
                    factor = particle->scalarDataSets[id].unitSI;

                    // Create the Y coordinates array
                    float *yarray = new float[particleBlock.numParticles];

                    // Read the dataset
                    openPMDFile.ReadParticleScalarBlock(yarray,
                                                        &factor,
                                                        H5T_FLOAT,
                                                        &particleBlock);

                    /* ____ Read the Z coordinates from the file _________ */

                    float *zarray = 0;
                    if (ndims>2)
                    {
                        // Dataset Id
                        id = particle->positionsId[0];

                        // Get the block properties for this dataSet
                        particle->GetBlockProperties(id, this->numTasks,
                                                     domain,
                                                     &particleBlock);

                        // Dataset path
                        strcpy(bufferDataSetName,particleBlock.dataSetPath);

                        // Multiplication factor
                        factor = particle->scalarDataSets[id].unitSI;

                        // Create the Z coordinates array
                        zarray = new float[particleBlock.numParticles];

                        // Read the dataset
                        openPMDFile.ReadParticleScalarBlock(zarray,
                                                            &factor,
                                                            H5T_FLOAT,
                                                            &particleBlock);

                    }

                    // Create the vtkPoints object and copy points into it.
                    vtkPoints *points = vtkPoints::New();
                    points->SetNumberOfPoints(particleBlock.numParticles);
                    float *pts = (float *) points->GetVoidPointer(0);
                    float *xc = xarray;
                    float *yc = yarray;
                    float *zc = zarray;
                    if(ndims == 3)
                    {
                        for(int i = 0; i < particleBlock.numParticles; ++i)
                        {
                            *pts++ = *xc++;
                            *pts++ = *yc++;
                            *pts++ = *zc++;
                        }
                    }
                    else if(ndims == 2)
                    {
                        for(int i = 0; i < particleBlock.numParticles; ++i)
                        {
                            *pts++ = *xc++;
                            *pts++ = *yc++;
                            *pts++ = 0.;
                        }
                    }

                    // Create a vtkUnstructuredGrid to contain the point cells.
                    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
                    ugrid->SetPoints(points);
                    points->Delete();
                    ugrid->Allocate(particleBlock.numParticles);
                    vtkIdType onevertex;
                    for(int i = 0; i < particleBlock.numParticles; ++i)
                    {
                        onevertex = i;
                        ugrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
                    }

                    // Delete temporary arrays.
                    delete [] xarray;
                    delete [] yarray;
                    delete [] zarray;
                    return ugrid;


                }
                // Double precision float
                else if (dataSize == 8)
                {
                    // Create the X coordinates array
                    double *xarray = new double[particleBlock.numParticles];

                    // Read the dataset
                    openPMDFile.ReadParticleScalarBlock(xarray,
                                                        &factor,
                                                        H5T_FLOAT,
                                                         &particleBlock);

                    /* ____ Read the Y coordinates from the file _________ */

                    // Dataset Id
                    id = particle->positionsId[1];

                    // Get the block properties for this dataSet
                    particle->GetBlockProperties(id, this->numTasks,
                                                 domain,
                                                 &particleBlock);

                    // Dataset path
                    strcpy(bufferDataSetName,particleBlock.dataSetPath);

                    // Multiplication factor
                    factor = particle->scalarDataSets[id].unitSI;

                    // Create the Y coordinates array
                    double *yarray = new double[particleBlock.numParticles];

                    // Read the dataset
                    openPMDFile.ReadParticleScalarBlock(yarray,
                                                        &factor,
                                                        H5T_FLOAT,
                                                        &particleBlock);

                    /* ____ Read the Z coordinates from the file _________ */

                    double *zarray = 0;
                    if (ndims>2)
                    {
                        // Dataset Id
                        id = particle->positionsId[0];

                        // Get the block properties for this dataSet
                        particle->GetBlockProperties(id, this->numTasks,
                                                     domain,
                                                     &particleBlock);

                        // Dataset path
                        strcpy(bufferDataSetName,particleBlock.dataSetPath);

                        // Multiplication factor
                        factor = particle->scalarDataSets[id].unitSI;

                        // Create the Z coordinates array
                        zarray = new double[particleBlock.numParticles];

                        // Read the dataset
                        openPMDFile.ReadParticleScalarBlock(zarray,
                                                            &factor,
                                                            H5T_FLOAT,
                                                            &particleBlock);

                    }

                    // Create the vtkPoints object and copy points into it.
                    vtkPoints *points = vtkPoints::New();
                    points->SetDataTypeToDouble();
                    points->SetNumberOfPoints(particleBlock.numParticles);
                    double *pts = (double *) points->GetVoidPointer(0);
                    double *xc = xarray;
                    double *yc = yarray;
                    double *zc = zarray;
                    if(ndims == 3)
                    {
                        for(int i = 0; i < particleBlock.numParticles; ++i)
                        {
                            *pts++ = *xc++;
                            *pts++ = *yc++;
                            *pts++ = *zc++;
                        }
                    }
                    else if(ndims == 2)
                    {
                        for(int i = 0; i < particleBlock.numParticles; ++i)
                        {
                            *pts++ = *xc++;
                            *pts++ = *yc++;
                            *pts++ = 0.;
                        }
                    }

                    // Create a vtkUnstructuredGrid to contain the point cells.
                    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
                    ugrid->SetPoints(points);
                    points->Delete();
                    ugrid->Allocate(particleBlock.numParticles);
                    vtkIdType onevertex;
                    for(int i = 0; i < particleBlock.numParticles; ++i)
                    {
                        onevertex = i;
                        ugrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
                    }

                    // Delete temporary arrays.
                    delete [] xarray;
                    delete [] yarray;
                    delete [] zarray;
                    return ugrid;

                }


            }
            /* ___ Sequential treatments _________________________________ */
            else
            {

                 // Dataset Id
                id = particle->positionsId[2];

                // Read the number of particles/elements from dataSet
                // It may be different from particle->numParticles
                nnodes = particle->scalarDataSets[id].numElements;

                // data size (float or double)
                dataSize = particle->scalarDataSets[id].dataSize;

                // Depending on the data size of this first dataSet,
                // we suppose that all data have the same type.

                // Simple precision float
                if (dataSize == 4)
                {

                    // Dataset path
                    strcpy(bufferDataSetName,particle->scalarDataSets[id].path);

                    // Multiplication factor
                    factor = particle->scalarDataSets[id].unitSI;

                    // Read the X coordinates from the file.
                    float *xarray = new float[nnodes];

                    // Read the dataset
                    openPMDFile.ReadScalarDataSet(xarray,nnodes,
                                                  &factor,
                                                  H5T_FLOAT,
                                                  bufferDataSetName);

                    // Read the Y coordinates from the file.
                    float *yarray = new float[nnodes];

                    // Dataset Id
                    id = particle->positionsId[1];

                    // Read the number of particles/elements from dataSet
                    // It may be different from particle->numParticles
                    nnodes = particle->scalarDataSets[id].numElements;

                    // Dataset path
                    strcpy(bufferDataSetName,particle->scalarDataSets[id].path);

                    // Multiplication factor
                    factor = particle->scalarDataSets[id].unitSI;

                    // Read the dataset
                    openPMDFile.ReadScalarDataSet(yarray,
                                                  nnodes,
                                                  &factor,
                                                  H5T_FLOAT,
                                                  bufferDataSetName);

                    float *zarray = 0;
                    if (ndims>2)
                    {
                        // Read the Z coordinates from the file.
                        zarray = new float[nnodes];

                        // Dataset Id
                        id = particle->positionsId[0];

                        // Read the number of particles/elements from dataSet
                        // It may be different from particle->numParticles
                        nnodes = particle->scalarDataSets[id].numElements;

                        // Dataset path
                        strcpy(bufferDataSetName,
                               particle->scalarDataSets[id].path);

                        // Multiplication factor
                        factor = particle->scalarDataSets[id].unitSI;

                        // Read the dataset
                        openPMDFile.ReadScalarDataSet(zarray,
                                                      nnodes,
                                                      &factor,
                                                      H5T_FLOAT,
                                                      bufferDataSetName);
                    }

                    // Create the vtkPoints object and copy points into it.
                    vtkPoints *points = vtkPoints::New();
                    points->SetNumberOfPoints(nnodes);
                    float *pts = (float *) points->GetVoidPointer(0);
                    float *xc = xarray;
                    float *yc = yarray;
                    float *zc = zarray;
                    if(ndims == 3)
                    {
                        for(int i = 0; i < nnodes; ++i)
                        {
                            *pts++ = *xc++;
                            *pts++ = *yc++;
                            *pts++ = *zc++;
                        }
                    }
                    else if(ndims == 2)
                    {
                        for(int i = 0; i < nnodes; ++i)
                        {
                            *pts++ = *xc++;
                            *pts++ = *yc++;
                            *pts++ = 0.;
                        }
                    }

                    // Create a vtkUnstructuredGrid to contain the point cells.
                    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
                    ugrid->SetPoints(points);
                    points->Delete();
                    ugrid->Allocate(nnodes);
                    vtkIdType onevertex;
                    for(int i = 0; i < nnodes; ++i)
                    {
                        onevertex = i;
                        ugrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
                    }

                    // Delete temporary arrays.
                    delete [] xarray;
                    delete [] yarray;
                    delete [] zarray;
                    return ugrid;

                }
                else if (dataSize == 8)
                {

                    // Dataset path
                    strcpy(bufferDataSetName,particle->scalarDataSets[id].path);

                    // Multiplication factor
                    factor = particle->scalarDataSets[id].unitSI;

                    // Read the X coordinates from the file.
                    double *xarray = new double[nnodes];

                    // Read the dataset
                    openPMDFile.ReadScalarDataSet(xarray,nnodes,
                                                  &factor,
                                                  H5T_FLOAT,
                                                  bufferDataSetName);

                    // Read the Y coordinates from the file.
                    double *yarray = new double[nnodes];

                    // Dataset Id
                    id = particle->positionsId[1];

                    // Read the number of particles/elements from dataSet
                    // It may be different from particle->numParticles
                    nnodes = particle->scalarDataSets[id].numElements;

                    // Dataset path
                    strcpy(bufferDataSetName,particle->scalarDataSets[id].path);

                    // Multiplication factor
                    factor = particle->scalarDataSets[id].unitSI;

                    // Read the dataset
                    openPMDFile.ReadScalarDataSet(yarray,
                                                  nnodes,
                                                  &factor,
                                                  H5T_FLOAT,
                                                  bufferDataSetName);

                    double *zarray = 0;
                    if (ndims>2)
                    {
                        // Read the Z coordinates from the file.
                        zarray = new double[nnodes];

                        // Dataset Id
                        id = particle->positionsId[0];

                        // Dataset path
                        strcpy(bufferDataSetName,
                               particle->scalarDataSets[id].path);

                        // Read the number of particles/elements from dataSet
                        // It may be different from particle->numParticles
                        nnodes = particle->scalarDataSets[id].numElements;

                        // Multiplication factor
                        factor = particle->scalarDataSets[id].unitSI;

                        // Read the dataset
                        openPMDFile.ReadScalarDataSet(zarray,
                                                      nnodes,
                                                      &factor,
                                                      H5T_FLOAT,
                                                      bufferDataSetName);

                        // Create the vtkPoints object and copy points into it.
                        vtkPoints *points = vtkPoints::New();
                        points->SetDataTypeToDouble();
                        points->SetNumberOfPoints(nnodes);
                        double *pts = (double *) points->GetVoidPointer(0);
                        double *xc = xarray;
                        double *yc = yarray;
                        double *zc = zarray;
                        if(ndims == 3)
                        {
                            for(int i = 0; i < nnodes; ++i)
                            {
                                *pts++ = *xc++;
                                *pts++ = *yc++;
                                *pts++ = *zc++;
                            }
                        }
                        else if(ndims == 2)
                        {
                            for(int i = 0; i < nnodes; ++i)
                            {
                                *pts++ = *xc++;
                                *pts++ = *yc++;
                                *pts++ = 0.;
                            }
                        }

                        // Create a vtkUnstructuredGrid to contain the point cells.
                        vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
                        ugrid->SetPoints(points);
                        points->Delete();
                        ugrid->Allocate(nnodes);
                        vtkIdType onevertex;
                        for(int i = 0; i < nnodes; ++i)
                        {
                            onevertex = i;
                            ugrid->InsertNextCell(VTK_VERTEX, 1, &onevertex);
                        }

                        // Delete temporary arrays.
                        delete [] xarray;
                        delete [] yarray;
                        delete [] zarray;
                        return ugrid;

                    }

                }
            }

        }

    }

    // No mesh name that we recognize.
    EXCEPTION1(InvalidVariableException, meshname);
    return 0;
}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
//  Modifications:
//      Mathieu Lobet, Dec 13 2016
//      I added the reading of single and double precision float datasets
//
//      Mathieu Lobet, Nov 25 2016
//      I added the parallel reading for the particles
//
//      Mathieu Lobet, Feb  8 2016
//      I added the reading of the double precision data
// ****************************************************************************

vtkDataArray *
avtOpenPMDFileFormat::GetVar(int timestate, int domain, const char *varname)
{
    debug5 << "avtOpenPMDFileFormat::GetVar("
         << "timestate=" << timestate << ","
         << "domain=" << domain << ","
         << "varname=" << varname << ")" <<endl;

    int     i,j,k,l,m;
    int     numValues;
    int     err;
    int     dims[3];
    string  buffer;
    char    xCompBuffer[128];
    char    yCompBuffer[128];
    hid_t   datasetId;
    hid_t   datasetType;
    hid_t   datasetSpace;
    int     ndims;
    float * array;
    float   factor;
    bool    varFound;

    // ________________________________________________________
    // FIELDS
    // This part get the scalar datasets from the fields

    // Shortcut pointer to the fields
    vector<PMDField> * fields = &(openPMDFile.iterations[timestate].fields);

    // Iteration over the fields
    for (std::vector<PMDField>::iterator field = fields->begin() ;
         field != fields->end(); ++field)
    {

        buffer = "Fields/";
        buffer += field->name;

        // Determine which scalar dataset to return.
        if (strcmp(varname, buffer.c_str()) == 0)
        {

            if (field->geometry == "cartesian")
            {

                // We treat the file in parallel by reading
                // the scalar dataset by block
                if (parallel)
                {

                    // Structure to store the block properties
                    fieldBlockStruct fieldBlock;

                    // We get the block properties
                    err = field->GetBlockProperties(this->numTasks,
                                                    domain,
                                                    &fieldBlock);

                    // Number of nodes
                    numValues = fieldBlock.nbTotalNodes;

                    // Factor for SI units
                    factor = field->unitSI;

                    // Float simple precision
                    if (field->dataSize == 4)
                    {
                        // Allocate the return vtkFloatArray object
                        vtkFloatArray * vtkArray = vtkFloatArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        float *data = (float *)vtkArray->GetVoidPointer(0);
                        // Reading of the dataset block
                        err = openPMDFile.ReadFieldScalarBlock(data,
                                                            &factor,
                                                            field->dataClass,
                                                            &fieldBlock);

                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }

                    }
                    // Float double precision
                    else if (field->dataSize == 8)
                    {
                        // Allocate the return vtkDoubleArray object
                        vtkDoubleArray * vtkArray = vtkDoubleArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        double *data = (double *)vtkArray->GetVoidPointer(0);
                        // Reading of the dataset block
                        err = openPMDFile.ReadFieldScalarBlock(data,
                                                            &factor,
                                                            field->dataClass,
                                                            &fieldBlock);

                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }

                    }
                    else
                    {
                        debug5 << " Error in avtOpenPMDFileFormat::GetVar" << endl;
                        debug5 << " The data size is not recognized." << endl;
                    }

                }
                // Sequential reading of the dataset
                else
                {
                    // Number of values
                    numValues = field->GetNumValues();
                    // Factor for SI units
                    factor = field->unitSI;

                    // Float simple precision
                    if (field->dataSize == 4)
                    {
                        // Allocate the return vtkFloatArray object
                        vtkFloatArray * vtkArray = vtkFloatArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        float *data = (float *)vtkArray->GetVoidPointer(0);
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(data,numValues,
                            &factor,field->dataClass,field->datasetPath);

                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }


                    }
                    // Float double precision
                    else if (field->dataSize == 8)
                    {
                        // Allocate the return vtkFloatArray object
                        vtkDoubleArray * vtkArray = vtkDoubleArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        double *data = (double *)vtkArray->GetVoidPointer(0);
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(data,numValues,
                            &factor,field->dataClass,field->datasetPath);

                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }

                    }
                    else
                    {
                        debug5 << " Error in avtOpenPMDFileFormat::GetVar" << endl;
                        debug5 << " The data size is not recognized." << endl;
                    }
                }

            }
            else if (field->geometry == "thetaMode")
            {

                // We treat the file in parallel by reading
                // the scalar dataset by block
                if (parallel)
                {

                    // Structure to store the block properties
                    fieldBlockStruct fieldBlock;

                    // We get the block properties
                    err = field->GetBlockProperties(this->numTasks,
                                                  domain,
                                                  &fieldBlock);

                    cerr << " Error in avtOpenPMDFileFormat::GetVar"
                         << endl;
                    cerr << " The parallel reading of the theta"
                         << "mode is not implemented"<< endl;
                }
                else
                {
                    // Dimensions
                    dims[0] = field->nbNodes[2]; // z direction
                    dims[1] = field->nbNodes[1]; // r direction
                    dims[2] = field->thetaNbNodes; // Theta direction

                    // Factor for SI units
                    factor = field->unitSI;

                    // Values to be read from the dataset
                    numValues = dims[1]*dims[0]*field->nbNodes[0];

                    // Simple precision
                    if (field->dataSize == 4)
                    {

                        // Allocate the dataset array with the different modes
                        float *dataSetArray = new float [numValues];

                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(dataSetArray,
                                                            numValues,
                                                            &factor,
                                                            H5T_FLOAT,
                                                            field->datasetPath);

                        // Total number of Values in the vtkarray
                        int numValuesVtkArray = dims[0]*dims[1]*dims[2];

                        // Allocate the return vtkFloatArray object. Note that
                        // you can use vtkFloatArray, vtkDoubleArray,
                        // vtkUnsignedCharArray, vtkIntArray, etc.
                        vtkFloatArray * vtkArray = vtkFloatArray::New();
                        vtkArray->SetNumberOfTuples(numValuesVtkArray);
                        float *data = (float *) vtkArray->GetVoidPointer(0);

                        // Generate the final data array with the modes
                        err = field->ComputeArrayThetaMode(dataSetArray,data);

                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }

                        // Delete temporary arrays
                        delete [] dataSetArray;

                    }
                    // Double precision
                    else if (field->dataSize == 8)
                    {

                        // Allocate the dataset array with the different modes
                        double *dataSetArray = new double [numValues];

                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(dataSetArray,
                                                            numValues,
                                                            &factor,
                                                            H5T_FLOAT,
                                                            field->datasetPath);

                        // Total number of Values in the vtkarray
                        numValues = dims[0]*dims[1]*dims[2];

                        // Allocate the return vtkFloatArray object. Note that
                        // you can use vtkFloatArray, vtkDoubleArray,
                        // vtkUnsignedCharArray, vtkIntArray, etc.
                        vtkDoubleArray * vtkArray = vtkDoubleArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        double *data = (double *)vtkArray->GetVoidPointer(0);

                        // Generate the final data array with the modes
                        err = field->ComputeArrayThetaMode(dataSetArray,data);

                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }

                        // Delete temporary arrays
                        delete [] dataSetArray;
                    }
                    else
                    {
                        cerr << " Error in avtOpenPMDFileFormat::GetVar" << endl;
                        cerr << " The data size is not recognized." << endl;
                    }
                }
            }

        }
    }

    // _______________________________________________
    // FIELD group



    // Temporary pointer to the field object
    PMDField * fieldTmp;

    // Shortcut pointer to the field groups
    vector<fieldGroupStruct> * fieldGroups
                             = &(openPMDFile.iterations[timestate].fieldGroups);

    // Iteration over the field group for additional components
    // computed with multiple datasets
    for (std::vector<fieldGroupStruct>::iterator
         fieldGroup = fieldGroups->begin();
         fieldGroup != fieldGroups->end();
         ++fieldGroup)
    {

      strcpy(xCompBuffer,"Fields/");
      strcat(xCompBuffer,fieldGroup->name);
      strcat(xCompBuffer,"/x");

      strcpy(yCompBuffer,"Fields/");
      strcat(yCompBuffer,fieldGroup->name);
      strcat(yCompBuffer,"/y");

      // Determine which additional field component to return.
      // If we want the x or y component in theta mode
      if  ((strcmp(varname, xCompBuffer) == 0) ||
           (strcmp(varname, yCompBuffer) == 0))
      {

          // And we are in theta mode
          if (fieldGroup->geometry == "thetaMode")
          {
              // check if components r and theta are available
              // If True, we compute x and y from r and theta.
              if ((fieldGroup->thetaComponents[0]>-1) &&
                  (fieldGroup->thetaComponents[1]>-1))
              {

                  // Temporary field pointer to the field object corresponding
                  //  to the r component in the list of fields
                  fieldTmp = &(openPMDFile.iterations[timestate].
                        fields[fieldGroup->thetaComponents[0]]);

                  // Dimensions
                  dims[0] = fieldTmp->nbNodes[2]; // z direction
                  dims[1] = fieldTmp->nbNodes[1]; // r direction
                  dims[2] = fieldTmp->thetaNbNodes; // Theta direction

                  // Factor for SI units
                  factor = fieldTmp->unitSI;

                  // Values to be read from the dataset
                  numValues = dims[1]*dims[0]*fieldTmp->nbNodes[0];

                  // Total number of Values in the final vtkarray
                  int numValuesFinalArray = dims[0]*dims[1]*dims[2];

                  // Simple precision
                  if (fieldTmp->dataSize == 4)
                  {
                      // Allocate the dataset array with the different modes
                      float *dataSetArray = new float [numValues];

                      // Reading of the dataset
                      err = openPMDFile.ReadScalarDataSet(dataSetArray,
                                                        numValues,
                                                        &factor,
                                                        H5T_FLOAT,
                                                        fieldTmp->datasetPath);

                     // Allocate array for the r component
                     float *rDataArray = new float [numValuesFinalArray];

                     // Generate the r data array with the modes
                     err = fieldTmp->ComputeArrayThetaMode(dataSetArray,rDataArray);

                     // The fieldTmp pointer now points to the new theta field component
                     fieldTmp = &(openPMDFile.iterations[timestate].
                             fields[fieldGroup->thetaComponents[1]]);

                     // Factor for SI units
                     factor = fieldTmp->unitSI;

                     // Reading of the dataset for the theta field component
                     err = openPMDFile.ReadScalarDataSet(dataSetArray,
                                                       numValues,
                                                       &factor,
                                                       H5T_FLOAT,
                                                       fieldTmp->datasetPath);

                     // Allocate the return vtkFloatArray object. Note that
                     // you can use vtkFloatArray, vtkDoubleArray,
                     // vtkUnsignedCharArray, vtkIntArray, etc.
                     vtkFloatArray * vtkArray = vtkFloatArray::New();
                     vtkArray->SetNumberOfTuples(numValuesFinalArray);
                     float *xDataArray = (float *)vtkArray->GetVoidPointer(0);

                     // Generate the theta field data array with the modes
                     // We put this array in xDataArray to save some memory
                     err = fieldTmp->ComputeArrayThetaMode(dataSetArray,xDataArray);

                     // Creation of dtheta, discretization of theta
                     float theta = 0;
                     float dtheta = 4.*acos(0.)/(fieldTmp->thetaNbNodes-1);

		     // Computation of the x component
		     if (strcmp(varname, xCompBuffer) == 0)
		     {
                         for(k = 0; k < fieldTmp->thetaNbNodes; ++k) // Loop theta
			 {
			   // theta
			   theta = k*dtheta;
			   for(j = 0; j < fieldTmp->nbNodes[1]; ++j) // Loop r
			     {
			       for(i = 0; i < fieldTmp->nbNodes[2]; ++i) // Loop z
				 {
				   // Absolute indexes
				   l = i + (j + k*fieldTmp->nbNodes[1])*fieldTmp->nbNodes[2];

				   // Update of data
				   xDataArray[l] = cos(theta)*rDataArray[l]
				     - sin(theta)*xDataArray[l];
				 }
			     }
			 }
		     }
		     // computation of the y component
		     else if (strcmp(varname, yCompBuffer) == 0)
		     {
                         for(k = 0; k < fieldTmp->thetaNbNodes; ++k) // Loop theta
			 {
			   // Computation of theta
			   theta = k*dtheta;
			   for(j = 0; j < fieldTmp->nbNodes[1]; ++j) // Loop r
			   {
			     for(i = 0; i < fieldTmp->nbNodes[2]; ++i) // Loop z
			     {
				 // Absolute indexes
				 l = i + (j + k*fieldTmp->nbNodes[1])*fieldTmp->nbNodes[2];

				 // Update of data
				 xDataArray[l] = sin(theta)*rDataArray[l]
				   + cos(theta)*xDataArray[l];
			     }
			   }
			 }
		     }

                     // If no error, we return the array
                     if (err>=0)
                     {
                         return vtkArray;
                     }

                  }
                  // Double precision
                  else if (fieldTmp->dataSize == 8)
                  {

                      // Allocate the dataset array with the different modes
                      double *dataSetArray = new double [numValues];

                      // Reading of the dataset
                      err = openPMDFile.ReadScalarDataSet(dataSetArray,
                                                        numValues,
                                                        &factor,
                                                        H5T_FLOAT,
                                                        fieldTmp->datasetPath);

                     // Allocate array for the r component
                     double *rDataArray = new double [numValuesFinalArray];

                     // Generate the r data array with the modes
                     err = fieldTmp->ComputeArrayThetaMode(dataSetArray,rDataArray);

		     // The fieldTmp pointer now points to the new theta field component
                     fieldTmp = &(openPMDFile.iterations[timestate].
                             fields[fieldGroup->thetaComponents[1]]);

                     // Factor for SI units
                     factor = fieldTmp->unitSI;

                     // Reading of the dataset for the theta field component
                     err = openPMDFile.ReadScalarDataSet(dataSetArray,
                                                       numValues,
                                                       &factor,
                                                       H5T_FLOAT,
                                                       fieldTmp->datasetPath);

                     // Allocate the return vtkFloatArray object. Note that
                     // you can use vtkFloatArray, vtkDoubleArray,
                     // vtkUnsignedCharArray, vtkIntArray, etc.
                     vtkDoubleArray * vtkArray = vtkDoubleArray::New();
                     vtkArray->SetNumberOfTuples(numValuesFinalArray);
                     double *xDataArray = (double *)vtkArray->GetVoidPointer(0);

                     // Generate the theta field data array with the modes
                     // We put this array in xDataArray to save some memory
                     err = fieldTmp->ComputeArrayThetaMode(dataSetArray,xDataArray);

                     // Creation of dtheta, discretization of theta
                     double theta = 0;
                     double dtheta = 4.*acos(0.)/(fieldTmp->thetaNbNodes-1);

		     // Computation of the x component
		     if (strcmp(varname, xCompBuffer) == 0)
		     {
                         for(k = 0; k < fieldTmp->thetaNbNodes; ++k) // Loop theta
			 {
			   // theta
			   theta = k*dtheta;
			   for(j = 0; j < fieldTmp->nbNodes[1]; ++j) // Loop r
			     {
			       for(i = 0; i < fieldTmp->nbNodes[2]; ++i) // Loop z
				 {
				   // Absolute indexes
				   l = i + (j + k*fieldTmp->nbNodes[1])*fieldTmp->nbNodes[2];

				   // Update of data
				   xDataArray[l] = cos(theta)*rDataArray[l]
				     - sin(theta)*xDataArray[l];
				 }
			     }
			 }
		     }
		     // computation of the y component
		     else if (strcmp(varname, yCompBuffer) == 0)
		     {
                         for(k = 0; k < fieldTmp->thetaNbNodes; ++k) // Loop theta
			 {
			   // Computation of theta
			   theta = k*dtheta;
			   for(j = 0; j < fieldTmp->nbNodes[1]; ++j) // Loop r
			   {
			     for(i = 0; i < fieldTmp->nbNodes[2]; ++i) // Loop z
			     {
				 // Absolute indexes
				 l = i + (j + k*fieldTmp->nbNodes[1])*fieldTmp->nbNodes[2];

				 // Update of data
				 xDataArray[l] = sin(theta)*rDataArray[l]
				   + cos(theta)*xDataArray[l];
			     }
			   }
			 }
		     }

                     // If no error, we return the array
                     if (err>=0)
                     {
                         return vtkArray;
                     }

                  }
              }
          }
      }
    }



    // ________________________________________________________
    // PARTICLES
    // This part gets the scalar datasets from the particles

    // Shortcut pointer to the particles
    vector<PMDParticle> * particles =
                            &(openPMDFile.iterations[timestate].particles);

    // Iteration over the particles
    for (std::vector<PMDParticle>::iterator particle = particles->begin() ;
         particle != particles->end(); ++particle)
    {
        // Iteration over the scalar datasets
        for (i=0;i<particle->scalarDataSets.size();i++)
        {

            buffer = "Particles/";
            buffer += particle->scalarDataSets[i].name;

            // Determine which scalar dataset to return.
            if (strcmp(varname, buffer.c_str()) == 0)
            {

                /* ___ Parallel ___________________________________ */

                // We treat the file in parallel by reading the scalar
                // dataset by block
                if (parallel)
                {

                    // Structure to store the block properties
                    particleBlockStruct particleBlock;

                    // We get the block properties
                    err = particle->GetBlockProperties(i, numTasks, domain ,
                                                       &particleBlock);

                    // multiplication factor for SI units
                    factor=particle->scalarDataSets[i].unitSI;

                    // Allocate the return vtkFloatArray object. Note that
                    // you can use vtkFloatArray, vtkDoubleArray,
                    // vtkUnsignedCharArray, vtkIntArray, etc.
                    if (particle->scalarDataSets[i].dataSize==4)
                    {
                        vtkFloatArray * vtkArray = vtkFloatArray::New();
                        vtkArray->SetNumberOfTuples(
                                  particleBlock.numParticles);

                        float *data = (float *)vtkArray->GetVoidPointer(0);
                        // Read particles
                        err = openPMDFile.ReadParticleScalarBlock(data,
                            &factor,
                            particle->scalarDataSets[i].dataClass,
                            &particleBlock);
                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }
                    }
                    // Double precision
                    else if (particle->scalarDataSets[i].dataSize==8)
                    {
                        vtkDoubleArray * vtkArray = vtkDoubleArray::New();
                        vtkArray->SetNumberOfTuples(
                                  particleBlock.numParticles);
                        double *data = (double *)vtkArray->GetVoidPointer(0);
                        // Read particles
                        err = openPMDFile.ReadParticleScalarBlock(data,
                            &factor,
                            particle->scalarDataSets[i].dataClass,
                            &particleBlock);
                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }
                    }

                }
                /* ___ Sequential ___________________________________ */
                else
                {
                    // Number of particles
                    numValues = particle->scalarDataSets[i].numElements;
                    // multiplication factor for SI units
                    factor=particle->scalarDataSets[i].unitSI;

                    // Allocate the return vtkFloatArray object. Note that
                    // you can use vtkFloatArray, vtkDoubleArray,
                    // vtkUnsignedCharArray, vtkIntArray, etc.
                    if (particle->scalarDataSets[i].dataSize==4)
                    {
                        vtkFloatArray * vtkArray = vtkFloatArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        float *data = (float *)vtkArray->GetVoidPointer(0);
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(data,
                            numValues,
                            &factor,
                            particle->scalarDataSets[i].dataClass,
                            particle->scalarDataSets[i].path);
                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }
                    }
                    // Double precision
                    else if (particle->scalarDataSets[i].dataSize==8)
                    {
                        vtkDoubleArray * vtkArray = vtkDoubleArray::New();
                        vtkArray->SetNumberOfTuples(numValues);
                        double *data = (double *)vtkArray->GetVoidPointer(0);
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(data,
                            numValues,
                            &factor,
                            particle->scalarDataSets[i].dataClass,
                            particle->scalarDataSets[i].path);
                        // If no error, we return the array
                        if (err>=0)
                        {
                            return vtkArray;
                        }
                    }
                }
            }
        }

    }
    //
    // If you have a file format where variables don't apply (for example a
    // strictly polygonal format like the STL (Stereo Lithography) format,
    // then uncomment the code below.

    EXCEPTION1(InvalidVariableException, varname);

    //
    // If you do have a scalar variable, here is some code that may be helpful.
    //
    // int ntuples = XXX; // this is the number of entries in the variable.
    // vtkFloatArray *rv = vtkFloatArray::New();
    // rv->SetNumberOfTuples(ntuples);
    // for (int i = 0 ; i < ntuples ; i++)
    // {
    //      rv->SetTuple1(i, VAL);  // you must determine value for ith entry.
    // }
    //
    // return rv;
    //

    return 0;

}


// ****************************************************************************
//  Method: avtOpenPMDFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      timestate  The index of the timestate.  If GetNTimesteps returned
//                 'N' time steps, this is guaranteed to be between 0 and N-1.
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: mlobet -- generated by xml2avt
//  Creation:   Tue Oct 18 11:49:37 PDT 2016
//
//  Modifications:
//  - Nov 25 2016 - Mathieu Lobet - Add parallel reading for the particles
//  - Feb  8 2017 - Mathieu Lobet - Add double precision for the particles
//
// ****************************************************************************

vtkDataArray *
avtOpenPMDFileFormat::GetVectorVar(int timestate,
                                   int domain,
                                   const char *varname)
{
    debug5 << "avtOpenPMDFileFormat::GetVectorVar("
         << "timestate=" << timestate << ","
         << "varname=" << varname
         << ")" << endl;

    int     i;
    int     err;
    int     dataSize;
    int     scalarDataSetId;
    int     numValues;
    int     numComponents;
    string  buffer;
    float   factor;

    // ________________________________________________________
    // PARTICLES
    // This part get the scalar datasets from the particles

    // Shortcut pointer to the particles
    vector<PMDParticle> * particles =
                            &(openPMDFile.iterations[timestate].particles);

    // Iteration over the particles
    for (std::vector<PMDParticle>::iterator particle = particles->begin() ;
         particle != particles->end(); ++particle)
    {
        // Iteration over the vector datasets
        for (i=0;i<particle->vectorDataSets.size();i++)
        {

            // Current dataset name to be compared with varname
            buffer = "Particles/";
            buffer += particle->vectorDataSets[i].name;

            // If we find the correct dataset of name varname...
            if (strcmp(varname, buffer.c_str()) == 0)
            {

                // Number of components (2D or 3D)
                numComponents = particle->numDimsMomenta;

                /* ___ Parallel ___________________________________ */

                // We treat the file in parallel by reading the scalar
                // dataset by block
                if (parallel)
                {

                    // Structure to store the block properties
                    particleBlockStruct particleBlock;

                    // Reading of the first dataset
                    scalarDataSetId = particle->vectorDataSets[i].dataSetId[2];

                    // dataSize
                    dataSize = particle->scalarDataSets[scalarDataSetId].dataSize;

                    // Based on th first scalar data size,
                    // we suppose the same fo the others

                    // Simple precision
                    if (dataSize == 4)
                    {
                        // Read component 1 from the file.
                        float *comp1 = new float[particleBlock.numParticles];

                        // We get the block properties
                        err = particle->GetBlockProperties(scalarDataSetId,
                                                           this->numTasks,
                                                           domain,
                                                           &particleBlock);
                        // Multiplication factor
                        factor = particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Copy of the dataset path
                        strcmp(particleBlock.dataSetPath,
                               particle->scalarDataSets[scalarDataSetId].path);
                        // Reading of the dataset
                        err = openPMDFile.ReadParticleScalarBlock(comp1,
                                                                  &factor,
                                                                  H5T_FLOAT,
                                                                  &particleBlock);

                        // Read component 2 from the file.
                        float *comp2 = new float[particleBlock.numParticles];
                        // Reading of the first dataset
                        scalarDataSetId = particle->vectorDataSets[i].dataSetId[1];
                        // We get the block properties
                        err = particle->GetBlockProperties(scalarDataSetId,
                                                           this->numTasks,
                                                           domain,
                                                           &particleBlock);
                        // Multiplication factor
                        factor = particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Copy of the dataset path
                        strcmp(particleBlock.dataSetPath,
                               particle->scalarDataSets[scalarDataSetId].path);
                        // Reading of the dataset
                        err = openPMDFile.ReadParticleScalarBlock(comp2,
                                                                  &factor,
                                                                  H5T_FLOAT,
                                                                  &particleBlock);


                        // Read component 3 from the file.
                        float *comp3;
                        if(numComponents > 2)
                        {
                            // Read component 3 from the file.
                            comp3 = new float[particleBlock.numParticles];
                            // Reading of the first dataset
                            scalarDataSetId =
                                         particle->vectorDataSets[i].dataSetId[0];
                            // We get the block properties
                            err = particle->GetBlockProperties(scalarDataSetId,
                                                               this->numTasks,
                                                               domain,
                                                               &particleBlock);
                            // Multiplication factor
                            factor =
                                 particle->scalarDataSets[scalarDataSetId].unitSI;
                            // Copy of the dataset path
                            strcmp(particleBlock.dataSetPath,
                                  particle->scalarDataSets[scalarDataSetId].path);
                            // Reading of the dataset
                            err = openPMDFile.ReadParticleScalarBlock(comp3,
                                                                      &factor,
                                                                      H5T_FLOAT,
                                                                &particleBlock);
                        }

                        // Allocate the return vtkFloatArray object. Note that
                        // you can use vtkFloatArray, vtkDoubleArray,
                        // vtkUnsignedCharArray, vtkIntArray, etc.
                        vtkFloatArray *arr = vtkFloatArray::New();
                        arr->SetNumberOfComponents(3);
                        arr->SetNumberOfTuples(particleBlock.numParticles);
                        float *data = (float *)arr->GetVoidPointer(0);
                        float *c1 = comp1;
                        float *c2 = comp2;
                        float *c3 = comp3;
                        if(numComponents == 3)
                        {
                            for(int i = 0; i < particleBlock.numParticles; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = *c3++;
                            }
                        }
                        else if(numComponents == 2)
                        {
                            for(int i = 0; i < particleBlock.numParticles; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = 0.;
                            }
                        }
                        else
                        {
                            delete [] comp1;
                            delete [] comp2;
                            delete [] comp3;
                            arr->Delete();
                            EXCEPTION1(InvalidVariableException, varname);
                        }
                        // Delete temporary arrays.
                        delete [] comp1;
                        delete [] comp2;
                        delete [] comp3;
                        return arr;
                    }
                    // Double precision
                    else if (dataSize==8)
                    {
                        // Read component 1 from the file.
                        double *comp1 = new double[particleBlock.numParticles];

                        // We get the block properties
                        err = particle->GetBlockProperties(scalarDataSetId,
                                                           this->numTasks,
                                                           domain,
                                                           &particleBlock);
                        // Multiplication factor
                        factor = particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Copy of the dataset path
                        strcmp(particleBlock.dataSetPath,
                               particle->scalarDataSets[scalarDataSetId].path);
                        // Reading of the dataset
                        err = openPMDFile.ReadParticleScalarBlock(comp1,
                                                                  &factor,
                                                                  H5T_FLOAT,
                                                              &particleBlock);

                        // Read component 2 from the file.
                        double *comp2 = new double[particleBlock.numParticles];
                        // Reading of the first dataset
                        scalarDataSetId =
                                       particle->vectorDataSets[i].dataSetId[1];
                        // We get the block properties
                        err = particle->GetBlockProperties(scalarDataSetId,
                                                           this->numTasks,
                                                           domain,
                                                           &particleBlock);
                        // Multiplication factor
                        factor =
                               particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Copy of the dataset path
                        strcmp(particleBlock.dataSetPath,
                               particle->scalarDataSets[scalarDataSetId].path);
                        // Reading of the dataset
                        err = openPMDFile.ReadParticleScalarBlock(comp2,
                                                                  &factor,
                                                                  H5T_FLOAT,
                                                              &particleBlock);

                        // Read component 3 from the file.
                        double *comp3;
                        if(numComponents > 2)
                        {
                            // Read component 3 from the file.
                            comp3 = new double[particleBlock.numParticles];
                            // Reading of the first dataset
                            scalarDataSetId =
                                         particle->vectorDataSets[i].dataSetId[0];
                            // We get the block properties
                            err = particle->GetBlockProperties(scalarDataSetId,
                                                               this->numTasks,
                                                               domain,
                                                               &particleBlock);
                            // Multiplication factor
                            factor =
                              particle->scalarDataSets[scalarDataSetId].unitSI;
                            // Copy of the dataset path
                            strcmp(particleBlock.dataSetPath,
                                particle->scalarDataSets[scalarDataSetId].path);
                            // Reading of the dataset
                            err = openPMDFile.ReadParticleScalarBlock(comp3,
                                                                      &factor,
                                                                      H5T_FLOAT,
                                                                &particleBlock);
                        }

                        // Allocate the return vtkDoubleArray object. Note that
                        // you can use vtkFloatArray, vtkDoubleArray,
                        // vtkUnsignedCharArray, vtkIntArray, etc.
                        vtkDoubleArray *arr = vtkDoubleArray::New();
                        arr->SetNumberOfComponents(3);
                        arr->SetNumberOfTuples(particleBlock.numParticles);
                        double *data = (double *)arr->GetVoidPointer(0);
                        double *c1 = comp1;
                        double *c2 = comp2;
                        double *c3 = comp3;
                        if(numComponents == 3)
                        {
                            for(int i = 0; i < particleBlock.numParticles; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = *c3++;
                            }
                        }
                        else if(numComponents == 2)
                        {
                            for(int i = 0; i < particleBlock.numParticles; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = 0.;
                            }
                        }
                        else
                        {
                            delete [] comp1;
                            delete [] comp2;
                            delete [] comp3;
                            arr->Delete();
                            EXCEPTION1(InvalidVariableException, varname);
                        }
                        // Delete temporary arrays.
                        delete [] comp1;
                        delete [] comp2;
                        delete [] comp3;
                        return arr;

                    }
                    // Else data size is not recognized
                    else
                    {
                        cerr << " Error in avtOpenPMDFileFormat::GetVectorVar"
                             << endl;
                        cerr << " DataSize is not recognized"
                             << endl;
                    }

                }
                /* ___ Sequential ___________________________________ */
                else
                {

                    // Reading of the first dataset
                    scalarDataSetId =
                                  particle->vectorDataSets[i].dataSetId[2];

                    // Number of particles for this dataset
                    numValues =
                        particle->scalarDataSets[scalarDataSetId].numElements;

                    // data size, we suppose that this is the
                    // same for the other components
                    dataSize =
                             particle->scalarDataSets[scalarDataSetId].dataSize;

                    // Simple precision
                    if (dataSize == 4)
                    {

                        // Read component 1 from the file.
                        float *comp1 = new float[numValues];

                        // Multiplication factor
                        factor =
                               particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(comp1,
                                                            numValues,
                                                            &factor,
                                                            H5T_FLOAT,
                                particle->scalarDataSets[scalarDataSetId].path);

                        // Read component 2 from the file.
                        float *comp2 = new float[numValues];
                        // Reading of the first dataset
                        scalarDataSetId =
                                      particle->vectorDataSets[i].dataSetId[1];
                        // Number of particles for this dataset
                        numValues =
                            particle->scalarDataSets[scalarDataSetId].numElements;

                        // Multiplication factor
                        factor =
                               particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(comp2,
                                                            numValues,
                                                            &factor,
                                                            H5T_FLOAT,
                              particle->scalarDataSets[scalarDataSetId].path);

                        float *comp3;
                        if(numComponents > 2)
                        {
                            // Read component 3 from the file.
                            comp3 = new float[numValues];
                            // Reading of the first dataset
                            scalarDataSetId =
                                      particle->vectorDataSets[i].dataSetId[0];
                            // Number of particles for this dataset
                            numValues =
                          particle->scalarDataSets[scalarDataSetId].numElements;
                            // Multiplication factor
                            factor =
                              particle->scalarDataSets[scalarDataSetId].unitSI;
                            // Reading of the dataset
                            err = openPMDFile.ReadScalarDataSet(comp3,
                                                                numValues,
                                                                &factor,
                                                                H5T_FLOAT,
                              particle->scalarDataSets[scalarDataSetId].path);
                        }

                        // Allocate the return vtkFloatArray object. Note that
                        // you can use vtkFloatArray, vtkDoubleArray,
                        // vtkUnsignedCharArray, vtkIntArray, etc.
                        vtkFloatArray *arr = vtkFloatArray::New();
                        arr->SetNumberOfComponents(3);
                        arr->SetNumberOfTuples(numValues);
                        float *data = (float *)arr->GetVoidPointer(0);
                        float *c1 = comp1;
                        float *c2 = comp2;
                        float *c3 = comp3;
                        if(numComponents == 3)
                        {
                            for(int i = 0; i < numValues; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = *c3++;
                            }
                        }
                        else if(numComponents == 2)
                        {
                            for(int i = 0; i < numValues; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = 0.;
                            }
                        }
                        else
                        {
                            delete [] comp1;
                            delete [] comp2;
                            delete [] comp3;
                            arr->Delete();
                            EXCEPTION1(InvalidVariableException, varname);
                        }
                        // Delete temporary arrays.
                        delete [] comp1;
                        delete [] comp2;
                        delete [] comp3;
                        return arr;
                    }
                    // Double precision
                    else if (dataSize==8)
                    {
                        // Read component 1 from the file.
                        double *comp1 = new double[numValues];
                        // Reading of the first dataset
                        scalarDataSetId =
                                      particle->vectorDataSets[i].dataSetId[2];
                        // Multiplication factor
                        factor =
                               particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(comp1,
                                                            numValues,
                                                            &factor,
                                                            H5T_FLOAT,
                                particle->scalarDataSets[scalarDataSetId].path);

                        // Read component 2 from the file.
                        double *comp2 = new double[numValues];
                        // Reading of the first dataset
                        scalarDataSetId =
                                      particle->vectorDataSets[i].dataSetId[1];
                        // Number of particles for this dataset
                        numValues =
                          particle->scalarDataSets[scalarDataSetId].numElements;
                        // Multiplication factor
                        factor =
                               particle->scalarDataSets[scalarDataSetId].unitSI;
                        // Reading of the dataset
                        err = openPMDFile.ReadScalarDataSet(comp2,
                                                            numValues,
                                                            &factor,
                                                            H5T_FLOAT,
                              particle->scalarDataSets[scalarDataSetId].path);

                        double *comp3;
                        if(numComponents > 2)
                        {
                            // Read component 3 from the file.
                            comp3 = new double[numValues];
                            // Reading of the first dataset
                            scalarDataSetId =
                                      particle->vectorDataSets[i].dataSetId[0];
                            // Number of particles for this dataset
                            numValues =
                         particle->scalarDataSets[scalarDataSetId].numElements;
                            // Multiplication factor
                            factor =
                              particle->scalarDataSets[scalarDataSetId].unitSI;
                            // Reading of the dataset
                            err = openPMDFile.ReadScalarDataSet(comp3,
                                                                numValues,
                                                                &factor,
                                                                H5T_FLOAT,
                              particle->scalarDataSets[scalarDataSetId].path);
                        }

                        // Allocate the return vtkDoubleArray object. Note that
                        // you can use vtkFloatArray, vtkDoubleArray,
                        // vtkUnsignedCharArray, vtkIntArray, etc.
                        vtkDoubleArray *arr = vtkDoubleArray::New();
                        arr->SetNumberOfComponents(3);
                        arr->SetNumberOfTuples(numValues);
                        double *data = (double *)arr->GetVoidPointer(0);
                        double *c1 = comp1;
                        double *c2 = comp2;
                        double *c3 = comp3;
                        if(numComponents == 3)
                        {
                            for(int i = 0; i < numValues; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = *c3++;
                            }
                        }
                        else if(numComponents == 2)
                        {
                            for(int i = 0; i < numValues; ++i)
                            {
                                *data++ = *c1++;
                                *data++ = *c2++;
                                *data++ = 0.;
                            }
                        }
                        else
                        {
                            delete [] comp1;
                            delete [] comp2;
                            delete [] comp3;
                            arr->Delete();
                            EXCEPTION1(InvalidVariableException, varname);
                        }
                        // Delete temporary arrays.
                        delete [] comp1;
                        delete [] comp2;
                        delete [] comp3;
                        return arr;

                    }
                    // The data size is not recognized
                    else
                    {
                      cerr << " Error in avtOpenPMDFileFormat::GetVectorVar"
                           << endl;
                      cerr << " DataSize is not recognized"
                           << endl;
                    }
                }
            }
        }
    }


    //
    // If you have a file format where variables don't apply (for example a
    // strictly polygonal format like the STL (Stereo Lithography) format,
    // then uncomment the code below.

    EXCEPTION1(InvalidVariableException, varname);


    return 0;

    //
    // If you do have a vector variable, here is some code that may be helpful.
    //
    // int ncomps = YYY;  // This is the rank of the vector - typically 2 or 3.
    // int ntuples = XXX; // this is the number of entries in the variable.
    // vtkFloatArray *rv = vtkFloatArray::New();
    // int ucomps = (ncomps == 2 ? 3 : ncomps);
    // rv->SetNumberOfComponents(ucomps);
    // rv->SetNumberOfTuples(ntuples);
    // float *one_entry = new float[ucomps];
    // for (int i = 0 ; i < ntuples ; i++)
    // {
    //      int j;
    //      for (j = 0 ; j < ncomps ; j++)
    //           one_entry[j] = ...
    //      for (j = ncomps ; j < ucomps ; j++)
    //           one_entry[j] = 0.;
    //      rv->SetTuple(i, one_entry);
    // }
    //
    // delete [] one_entry;
    // return rv;
    //
}
